#ifndef MATRIX_HPP
#define MATRIX_HPP

namespace TermStructure {
    //! A matrix using double pointer as the storage.
	template <typename T, typename TSize = long>
	class SimpleMatrix
	{
        //! Sizes
		TSize _length1;
		TSize _length2;
        //! the storage of the matrix
		T** _storage;
        //! initate the matrix
		void _init(const T &defaultT = 0)
		{
			_storage = new T*[_length1];
			for(TSize i = 0; i < _length1; i++)
			{
				_storage[i] = new T[_length2];
				for(TSize j = 0; j < _length2; j++)
				{
					_storage[i][j] = defaultT;
				}
			}
		}
        //! Release all resources
		void _delocate()
		{
			for(TSize i = 0; i < _length1; i++)
			{
				delete[] _storage[i];
			}
			delete[] _storage;
		}
	public:
		//! Constructors
		SimpleMatrix(TSize length1 = 1, TSize length2 = 1, const T &defaultT = 0):_length1(length1),_length2(length2)
		{
			_init(defaultT);
		}
        //! copy constructor
		SimpleMatrix(const SimpleMatrix &anotherMatrix):_length1(anotherMatrix._length1), _length2(anotherMatrix._length2)
		{
			_init();
			for(TSize i = 0; i < _length1; i++)
			{
				_storage[i] = new T[_length2];
				for(TSize j = 0; j < _length2; j++)
				{
					_storage[i][j] = anotherMatrix(i,j);
				}
			}            
		}
		//! operator=
		SimpleMatrix &operator=(const SimpleMatrix &anotherMatrix)
		{
			if(&anotherMatrix==this) return *this;
			if(anotherMatrix._length1 != _length1 || anotherMatrix._length2 != _length2)
			{
				this->_delocate();
			}
			for(TSize i = 0; i < _length1; i++)
			{
				_storage[i] = new T[_length2];
				for(TSize j = 0; j < _length2; j++)
				{
					_storage[i][j] = anotherMatrix(i,j);
				}
			}
			return *this;
		}
        //! Resize the matrix. The matrix's contents will be released
		SimpleMatrix &Resize(const TSize &i, const TSize &j)
		{
			if(this->_length1 == i && this->_length2 == j)
			{
				return *this;
			}
			else
			{
				_delocate();
				this->_length1 = i;
				this->_length2 = j;
				_init();
			}
			return *this;
		}
        //! Get the length
		TSize GetLength(int i) const
		{
			if(i == 0) return _length1;
			else if(i==1) return _length2;
			return -1;
		}
        //! DEstructor
		virtual ~SimpleMatrix()
		{
			_delocate();
		}
        //! Get the i,j element
		T &operator()(TSize i, TSize j)
		{
			return _storage[i][j];
		}
        //! Get the i,j element
		const T &operator()(TSize i, TSize j) const
		{
			return _storage[i][j];
		}
		
	};
    //! Wrap a C++ matrix, T[][] to a matrix satisfies the requirements of the matrix
	template<typename T, int NPrimary, int NSecondary, typename TSize = long>
	class MatrixWrapper
	{
        //! underlying C++ matrix
		T (*_underlying)[NPrimary][NSecondary];
	public:
        //! Default Constructor
		MatrixWrapper(T (&underlying)[NPrimary][NSecondary]):_underlying(&underlying)
		{
		}
        //! Get the i,j element
		T &operator()(TSize i, TSize j)
		{
			return (*_underlying)[i][j];
		}
        //! Get the i,j element
		const T &operator()(TSize i, TSize j) const
		{
			return (*_underlying)[i][j];
		}
        //! Get the length
		TSize GetLength(int i) const
		{
			if(i==0) return NPrimary;
			else return NSecondary;
		}
	};
}



#endif
